function [A,b] = assemble_Elliptic(dofs,V,T,f,caseNum)
% shapes and shapes_deri
  d = max(T(:,4));
%   [a,c]=get_demo_case(caseNum);
  a=1;c=0;
  m = (d+1)*(d+2)/2 ;
  nt = size(T, 1);
  dim = max(max(dofs));
  b = zeros(dim, 1);
  A = spalloc(dim,dim,m*m*nt); 
  [qw,qp] = quad_rule(13);
  % reserve spaces
  shape_value_cell = cell(10,1);
  shape_deri_x_cell = cell(10,1);
  shape_deri_y_cell = cell(10,1);
  % pre-set the value of the cell matrix
  for i=1:10
      [shape_value_cell{i},shape_deri_x_cell{i},shape_deri_y_cell{i}] = eval_bb(i,qp(1,:),qp(2,:));
  end
  % note that must be evaluate at master element ,so use qp(:,:);
  for tri = 1:nt
      d = T(tri,4);
      m = (d+1)*(d+2)/2 ;
      shape_value = shape_value_cell{d};
      shape_deri_x = shape_deri_x_cell{d};
      shape_deri_y = shape_deri_y_cell{d};
      v1=V(T(tri,1),:);v2=V(T(tri,2),:);v3=V(T(tri,3),:);
      x13 = v1(1) - v3(1);  %a
      x23 = v2(1) - v3(1);  %-c
      y23 = v2(2) - v3(2);  %d
      y13 = v1(2) - v3(2);  %-b
      J = x13*y23 - x23*y13;  
      a1 = (x23^2 + y23^2)/J^2;  
      a12 = (x23*x13 + y23*y13)/J^2; 
      a2 = (x13^2 + y13^2)/J^2;  
  
      qp_glob = zeros(size(qp));
      qp_glob(1,:) = v3(1) + x13*qp(1,:) + x23*qp(2,:);
      qp_glob(2,:) = v3(2) + y13*qp(1,:) + y23*qp(2,:);  
      % evaluate at integrate points, 
      fvec = feval(f, qp_glob(1,:), qp_glob(2,:),caseNum);
  
  %  then assemble the system:
      A_loc = zeros(m,m);
      b_loc = zeros(m,1);    % b_loc = J*(((ones(m,1)*fvec).*shapes(qp))*(qw'));
      for i=1:m
         for j=1:m
            mass = shape_value(i,:).*shape_value(j,:)*c ;
            stiff = (a1*shape_deri_x(i,:).*shape_deri_x(j,:)- ...
                  a12*(shape_deri_x(i,:).*shape_deri_y(j,:)+shape_deri_y(i,:).*shape_deri_x(j,:))+...
                  a2*shape_deri_y(i,:).*shape_deri_y(j,:))*a;
            A_loc(i,j) = (mass+stiff)*(qw')*J;
         end
         b_loc(i) = (fvec.*shape_value(i,:))*(qw')*J;
      end  
      if(size(dofs,2)==2)
          loc_dof = dofs(tri,1):dofs(tri,2);
      else
          loc_dof = dofs(tri,:);
      end
      A(loc_dof,loc_dof) = A(loc_dof,loc_dof) + A_loc;
      b(loc_dof) =  b(loc_dof) + b_loc; 
  end
